
package Matrices;

import main.Errores;

/**
  *  4g3nt
  *  Copyright (C) 2008  Carlos Molinero y César Andrés
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
  *  the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU General Public License for more details.
  *
  *    You should have received a copy of the GNU General Public License
  *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  *
  */
public class Matriz {

    private double[][] _matriz;
    public Matriz(int dimension) 
    {
        _matriz=new double[dimension][dimension];
        for (int i=0; i < dimension; i++)
            for (int j=0; j < dimension; j++)
                _matriz[i][j]=0;
        }
    /**
     * Constructora por paso de doble vector
     * @param x
     */
    public Matriz(double[][] matriz) 
    {
        _matriz=matriz;
    }

    public Matriz (Matriz matriz)
    {
        int dimension = matriz.getDimension();
        _matriz=new double[dimension][dimension];
        for (int i=0; i < dimension; i++)
            for (int j=0; j < dimension; j++)
                _matriz[i][j]=matriz._matriz[i][j];
    }

    double traza()
    {
        double tr=0.0;
        for( int i = 0 ; i < getDimension() ; i++)
             tr+=_matriz[i][i];
        return tr;
    }
    /**
     * Suma de dos matrices cuadradas
     *  
     * Devuelve 
     * @param a
     * @param b
     * @return
     */
    public static Matriz suma(Matriz a, Matriz b)
    {
        Matriz resultado=new Matriz(a.getDimension());
        for(int i = 0 ; i < a.getDimension(); i++){
            for(int j = 0 ;  j < a.getDimension(); j++)
                resultado._matriz[i][j]=a._matriz[i][j]+b._matriz[i][j];
            
        }
        return resultado;
    }
    /**
     * Producto de dos matrices
     * @param a
     * @param b
     * @return
     */
    public static Matriz producto(Matriz a, Matriz b)
    {
        int n = a.getDimension();
        Matriz resultado=new Matriz(n);
        for(int i = 0 ; i <n ; i++){
            for(int j = 0 ; j <n ; j++)
                for(int k = 0 ; k < n; k++)
                    resultado._matriz[i][j]+=a._matriz[i][k]*b._matriz[k][j];
        }
        return resultado;
    }
    /**
     * Producto de una matriz por un escalar
     * @param a
     * @param d
     * @return
     */
    public static Matriz producto(Matriz a, double d)
    {
        int n = a.getDimension();
        Matriz resultado = new Matriz(n);
        for(int i = 0 ;  i< n ; i++)
            for(int j = 0 ; j <n ; j++)
                resultado._matriz[i][j]=a._matriz[i][j]*d;
         return resultado;
    }
    /**
     * Producto de un escalar por una matriz
     * @param d
     * @param a
     * @return
     */
    public static Matriz producto(double d, Matriz a)
    {
        return producto(a,d);
    }
    /**
     * Producto de un vector fila por una matriz da un vector fila 
     *          (1xn) (nxn)= (1xn)
     * @param v
     * @param a
     * @return
     */
    public static Vector producto(Vector v, Matriz a)
    {
        int dimension=v.getDimension();  //dimensión
        double[] b=new double[dimension];
        for (int i = 0 ; i < dimension ; i++)
            for (int j = 0; j < dimension; j++)
                b[i] = v.getPosicion(j)*a._matriz[j][i];
        
        return new Vector(b);
    }
    /**
     * producto de una matriz por un vector columna 
     *      (nxn) (nx1)= (nx1)
     * @param a
     * @param v
     * @return
     */
    public static Vector producto(Matriz a, Vector v)
    {
        int n=v.getDimension();  //dimensión
        double[] b=new double[n];
        for(int i = 0 ; i < n ;  i++)
            for(int j = 0 ; j < n ; j++)
                b[i]+=a._matriz[i][j]*v.getPosicion(j);
        return new Vector(b);
    }
    /**
     * Determinante de la matriz
     * @return
     */
    public double determinante()
    {
        
        int n = getDimension();
        Matriz a= new Matriz(this);
        for(int k = 0 ; k < n-1 ; k++)
            for(int i = k+1 ; i < n ; i++)
                for(int j=k+1; j<n; j++)
                    a._matriz[i][j]-=a._matriz[i][k]*a._matriz[k][j]/a._matriz[k][k];
        double deter=1.0;
        for(int i = 0 ; i < n; i++)  deter*=a._matriz[i][i];
        return deter;
    }
    /**
     * matriz inversa
     * @param d
     * @return
     */
    public static Matriz inversa(Matriz d){
        int n=d.getDimension();  //dimensión de la matriz
        Matriz a=new Matriz(d);
        Matriz b=new Matriz(n);   //matriz de los términos independientes
        Matriz c=new Matriz(n);   //matriz de las incógnitas
        //matriz unidad
        for(int i = 0 ; i < n ; i++){
            b._matriz[i][i]=1.0;
        }
        //transformación de la matriz y de los términos independientes
        for(int k = 0 ; k < n-1 ; k++)
            for(int i = k+1; i < n; i++){
            //términos independientes
                for(int s=0; s<n; s++)
                    b._matriz[i][s]-=a._matriz[i][k]*b._matriz[k][s]/a._matriz[k][k];
            //elementos de la matriz
                for(int j=k+1; j<n; j++)
                    a._matriz[i][j]-=a._matriz[i][k]*a._matriz[k][j]/a._matriz[k][k];
                
            }
        //cálculo de las incógnitas, elementos de la matriz inversa
        for(int s = 0 ; s < n ; s++){
            c._matriz[n-1][s]=b._matriz[n-1][s]/a._matriz[n-1][n-1];
            for(int i = n-2 ; i >= 0 ; i--){
                c._matriz[i][s]=b._matriz[i][s]/a._matriz[i][i];
                for(int k=n-1; k>i; k--){
                    c._matriz[i][s]-=a._matriz[i][k]*c._matriz[k][s]/a._matriz[i][i];
                }
            }
        }
        return c;
    }

    /**
     * Matriz traspuesta
     * @param a
     * @return
     */
    public static Matriz traspuesta(Matriz a){
        int n=a.getDimension();    //dimensión
        Matriz d=new Matriz(n);
        for(int i = 0 ; i < n; i++)
            for(int j = 0 ; j < n ; j++)
                d._matriz[i][j]=a._matriz[j][i];
        return d;
    }
    /**
     * Polinomio característico
     * @return
     */
    public double[] polCaracteristico(){
        int n = getDimension();
        Matriz pot=new Matriz(n);
        //matriz unidad
        for(int i = 0 ; i < n ; i++)
            pot._matriz[i][i]=1.0;
        double[] p=new double[n+1];
        double[] s=new double[n+1];
        for(int i = 1 ; i <=n ; i++)
            {
                pot=Matriz.producto(pot, this);
                s[i]=pot.traza();
            }
        p[0]=1.0;
        p[1]=-s[1];
        for(int i = 2 ; i<=n ; i++)
            {
                p[i]=-s[i]/i;
                for(int j=1 ; j<i ; j++)
                    p[i]-=s[i-j]*p[j]/i;
            
            }
        return p;
    }
    
    public Matriz valoresPropios(double[] valores, int maxIter){
        final double CERO=1e-8;
        int n = getDimension();
        double maximo, tolerancia, sumsq;
        double x, y, z, c, s;
        int contador=0;
        int i, j, k, l;
        Matriz a=new Matriz(this);     
        //matriz copia
        Matriz p=new Matriz(n);
        Matriz q=new Matriz(n);
        //matriz unidad
        for(i=0 ; i < n ; i++)
            {
                q._matriz[i][i]=1.0;
            }
        while(contador<maxIter)
            {
            k=0; l=1;
            maximo=Math.abs(a._matriz[k][1]);
            for(i=0 ; i <n-1 ; i++)
                for(j=i+1 ; j < n ; j++)
                    if(Math.abs(a._matriz[i][j])>maximo){
                        k=i;        
                        l=j;
                        maximo=Math.abs(a._matriz[i][j]);
                    }
                
            
            sumsq=0.0;
            for(i=0 ; i < n ; i++)
                sumsq+=a._matriz[i][i]*a._matriz[i][i];
            
            tolerancia=0.0001*Math.sqrt(sumsq)/n;
            if(maximo<tolerancia) break;
            //calcula la matriz ortogonal de p
            //inicialmente es la matriz unidad
            for(i=0 ; i < n ; i++)
                for(j=0 ; j < n ; j++)
                    if (i==j) p._matriz[i][i]=1.0;
                    else p._matriz[i][j]=0.0;
            
            
            y=a._matriz[k][k]-a._matriz[l][l];
            if(Math.abs(y)<CERO)
                {
                c=s=Math.sin(Math.PI/4);
                }
                else
                {
                    x=2*a._matriz[k][l];
                    z=Math.sqrt(x*x+y*y);
                    c=Math.sqrt((z+y)/(2*z));
                    s=signo(x/y)*Math.sqrt((z-y)/(2*z));
                }
            
            p._matriz[k][k]=c;
            p._matriz[l][l]=c;
            p._matriz[k][l]=s;
            p._matriz[l][k]=-s;
            a=Matriz.producto(p, Matriz.producto(a, Matriz.traspuesta(p)));
            q=Matriz.producto(q, Matriz.traspuesta(p));
            contador++;
        }

        if(contador==maxIter)
            Errores.insertar("No se han podido calcular los valores propios");
        
        //valores propios
        //double[] valores=new double[n];
        for(i=0; i < n ; i++)
            valores[i]=(double)Math.round(a._matriz[i][i]*1000)/1000;
        //vectores propios
        return q;
    }

    private int signo(double x)
        {
        return (x>0 ? 1 : -1);
        }

    @Override
    public String toString()
    {
        String texto="\n";
        int n = getDimension();
        for(int i=0 ;  i < n ; i++){
            for(int j = 0 ; j < n ; j++)
                texto+="\t "+(double)Math.round(1000*_matriz[i][j])/1000;
            texto+="\n";
        }
        texto+="\n";
        return texto;
  }
    public int getDimension()
    {
    return _matriz[0].length;
    }

}